function [ outEFG ] = wrapperEFGArray( xNodes, xSamples, optEFG)
%WRAPPEREFG Returns a wrapper of the EFG aproximants and associated 
%   gradients.
%
%   INPUT:
%   xNodes      :   Set of nodes.
%   xSamples    :   Set of samples points.
%   optEFG      :   Wrapper structure that contains rho, poliBaseType 
%                   and order fields.  
%
%   OUTPUT
%   outEFG.p_samp   :   Shape functions at each sample point.
%   outEFG.dp_samp  :   Gradient of the shape functions at each sample point.


rho             = optEFG.rho;
poliBaseType    = optEFG.poliBaseType;
order           = optEFG.order;   

totalNodes      = size(xNodes,1);
totalSamples    = size(xSamples,1);
dimension       = size(xSamples,2);
baseDimension   = getPolibaseDimension(order);

samplesShapeFunction = zeros(totalSamples,totalNodes);
gradSamplesShapeFunction = zeros(totalSamples,totalNodes,dimension);

parfor indexSample = 1 : totalSamples

    %   Inner loop initialization
    gradAlpha = zeros(baseDimension, dimension);
    gradN     = zeros(totalNodes, dimension);
    weights   = zeros(totalNodes, 1);

    for indexNode = 1 : totalNodes
        distance = norm(xSamples(indexSample, : ) - xNodes(indexNode, : ));
        weights(indexNode) = gaussianWeight(distance , rho);
    end
    B = assemblyBMatrix( xSamples(indexSample, : ), xNodes, order , weights , poliBaseType );
    A = assemblyMomentMatrix( xSamples(indexSample, : ), xNodes, weights , poliBaseType, order );
    p = poliBase(xSamples(indexSample, : ), order, poliBaseType);
    [alpha L U]= solveAlpha (A,p);
    N = alpha * B';

    samplesShapeFunction(indexSample, : ) = N;

    gradP       = gradPoliBase (xSamples(indexSample, : ), order, poliBaseType);
    gradWeights = gradGaussianWeight (xNodes,xSamples(indexSample, : ), rho);
    gradA       = assemblyGradMomentMatrix (xSamples(indexSample, : ), xNodes,...
                  gradWeights, poliBaseType, order);
    loadVectorGradAlpha = assemblyLoadVectorGradAlpha ...
                        (gradP(:,1), gradP (:,2),...
                        gradA(:,:,1), gradA(:,:,2), alpha );

    for d = 1 : dimension
        gradAlpha(:,d) = solveGradAlpha(L, U,loadVectorGradAlpha(:,d)');
    end

    gradB = assemblyGradBMatrix(xSamples(indexSample, : ), xNodes, order,...
        gradWeights, poliBaseType);

    for d = 1 : dimension
        gradN(:,d) = gradAlpha(:,d)' * B' +  alpha * gradB(:,:,d)';
    end
    
    gradSamplesShapeFunction(indexSample, : , : ) = gradN; 
end


outEFG.p_samp = samplesShapeFunction;
outEFG.dp_samp = gradSamplesShapeFunction;

end

